lavaanパッケージ

lavaanパッケージは、Rで共分散構造分析(structual equation modeling)を実行するパッケージです。公式サイトは以下の通りです:

公式サイトにある確認的因子分析とSEMを実行すると、以下のとおりとなる:

library(lavaan)
## This is lavaan 0.5-20
## lavaan is BETA software! Please report any bugs.
xx <- HolzingerSwineford1939

HS_model <- '
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
'

fit <- cfa(HS_model, data=xx)
summary(fit)
## lavaan (0.5-20) converged normally after  35 iterations
## 
##   Number of observations                           301
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic               85.306
##   Degrees of freedom                                24
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  Z-value  P(>|z|)
##   visual =~                                           
##     x1                1.000                           
##     x2                0.554    0.100    5.554    0.000
##     x3                0.729    0.109    6.685    0.000
##   textual =~                                          
##     x4                1.000                           
##     x5                1.113    0.065   17.014    0.000
##     x6                0.926    0.055   16.703    0.000
##   speed =~                                            
##     x7                1.000                           
##     x8                1.180    0.165    7.152    0.000
##     x9                1.082    0.151    7.155    0.000
## 
## Covariances:
##                    Estimate  Std.Err  Z-value  P(>|z|)
##   visual ~~                                           
##     textual           0.408    0.074    5.552    0.000
##     speed             0.262    0.056    4.660    0.000
##   textual ~~                                          
##     speed             0.173    0.049    3.518    0.000
## 
## Variances:
##                    Estimate  Std.Err  Z-value  P(>|z|)
##     x1                0.549    0.114    4.833    0.000
##     x2                1.134    0.102   11.146    0.000
##     x3                0.844    0.091    9.317    0.000
##     x4                0.371    0.048    7.779    0.000
##     x5                0.446    0.058    7.642    0.000
##     x6                0.356    0.043    8.277    0.000
##     x7                0.799    0.081    9.823    0.000
##     x8                0.488    0.074    6.573    0.000
##     x9                0.566    0.071    8.003    0.000
##     visual            0.809    0.145    5.564    0.000
##     textual           0.979    0.112    8.737    0.000
##     speed             0.384    0.086    4.451    0.000
model2 <- '
  # measurement model
    ind60 =~ x1 + x2 + x3
    dem60 =~ y1 + y2 + y3 + y4
    dem65 =~ y5 + y6 + y7 + y8
  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60
  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 + y6
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'
fit2 <- sem(model2, data=PoliticalDemocracy)
summary(fit2)
## lavaan (0.5-20) converged normally after  68 iterations
## 
##   Number of observations                            75
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic               38.125
##   Degrees of freedom                                35
##   P-value (Chi-square)                           0.329
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Standard Errors                             Standard
## 
## Latent Variables:
##                    Estimate  Std.Err  Z-value  P(>|z|)
##   ind60 =~                                            
##     x1                1.000                           
##     x2                2.180    0.139   15.742    0.000
##     x3                1.819    0.152   11.967    0.000
##   dem60 =~                                            
##     y1                1.000                           
##     y2                1.257    0.182    6.889    0.000
##     y3                1.058    0.151    6.987    0.000
##     y4                1.265    0.145    8.722    0.000
##   dem65 =~                                            
##     y5                1.000                           
##     y6                1.186    0.169    7.024    0.000
##     y7                1.280    0.160    8.002    0.000
##     y8                1.266    0.158    8.007    0.000
## 
## Regressions:
##                    Estimate  Std.Err  Z-value  P(>|z|)
##   dem60 ~                                             
##     ind60             1.483    0.399    3.715    0.000
##   dem65 ~                                             
##     ind60             0.572    0.221    2.586    0.010
##     dem60             0.837    0.098    8.514    0.000
## 
## Covariances:
##                    Estimate  Std.Err  Z-value  P(>|z|)
##   y1 ~~                                               
##     y5                0.624    0.358    1.741    0.082
##   y2 ~~                                               
##     y4                1.313    0.702    1.871    0.061
##     y6                2.153    0.734    2.934    0.003
##   y3 ~~                                               
##     y7                0.795    0.608    1.308    0.191
##   y4 ~~                                               
##     y8                0.348    0.442    0.787    0.431
##   y6 ~~                                               
##     y8                1.356    0.568    2.386    0.017
## 
## Variances:
##                    Estimate  Std.Err  Z-value  P(>|z|)
##     x1                0.082    0.019    4.184    0.000
##     x2                0.120    0.070    1.718    0.086
##     x3                0.467    0.090    5.177    0.000
##     y1                1.891    0.444    4.256    0.000
##     y2                7.373    1.374    5.366    0.000
##     y3                5.067    0.952    5.324    0.000
##     y4                3.148    0.739    4.261    0.000
##     y5                2.351    0.480    4.895    0.000
##     y6                4.954    0.914    5.419    0.000
##     y7                3.431    0.713    4.814    0.000
##     y8                3.254    0.695    4.685    0.000
##     ind60             0.448    0.087    5.173    0.000
##     dem60             3.956    0.921    4.295    0.000
##     dem65             0.172    0.215    0.803    0.422

SEMでの結果は、多くの場合ダイアグラムで図示することが多いけれど、このlavaanパッケージはplotを自動的には出力しません。

DiagrammeRパッケージ

Rでダイアグラムを描くパッケージです。公式サイトは以下のとおりです:

公式サイトにあるサンプルコードを実行すると、以下のとおりです:

library(DiagrammeR)
# Create a node data frame
nodes <-
  create_nodes(
    nodes = c("a", "b", "c", "d"),
    label = FALSE,
    type = "lower",
    style = "filled",
    color = "aqua",
    shape = c("circle", "circle",
              "rectangle", "rectangle"),
    data = c(3.5, 2.6, 9.4, 2.7)
  )

edges <-
  create_edges(
    from = c("a", "b", "c"),
    to = c("d", "c", "a"),
    rel = "leading_to"
  )

graph <-
  create_graph(
    nodes_df = nodes,
    edges_df = edges,
    node_attrs = "fontname = Helvetica",
    edge_attrs = c("color = blue",
                   "arrowsize = 2")
  )

render_graph(graph)

lavaan2DgmR開発の試み

この2つを利用し、lavaanの出力オブジェクトを自動的にDiagrammeRのgraphオブジェクトを生成する関数を作成しました。

create_graph_lavaan <- function(fit, ...) {
  library(Hmisc) #順番注意
  library(DiagrammeR)
  library(dplyr)
  
  # 必要データ取り出し
  d <- data.frame(fit@ParTable) %>% 
    select(lhs,op,rhs,est)
  #こっちもありかも
  #こっちを使うなら、推定値だけじゃなくてp値も持ってこれるので、
  #アスタリスクを自動で付けたりとかできるようになる
  #library(lavaan)
  #d <- fit %>% parameterEstimates() %>% 
  #  select(lhs,op,rhs,est)
  
  # make NDF
  latent <- d %>%
    filter(op == "=~") %>%
    select(nodes = lhs) %>%
    distinct %>%
    mutate(shape = "circle", rank = "same")
  observed <- d %>% 
    filter(op != "~1", lhs %nin% latent$nodes) %>%
    select(nodes = lhs) %>%
    distinct %>%
    mutate(shape = "square")
  nodes_df <- combine_nodes(latent, observed)
  # make EDF
  all_paths <- d %>%
    filter(op != "~1") %>%
    mutate(label = format(round(est, 2),nsmall=2)) %>%
    select(-est)
  loadings <- all_paths %>%
    filter(op == "=~") %>%
    mutate(from = lhs, to = rhs, style = "dashed") %>%
    select(from, to, style, label)
  regressions <- all_paths %>%
    filter(op == "~~") %>%
    mutate(to = lhs, from = rhs, style = "solid", dir = "both") %>%
    select(from, to, style, label, dir)
  predictions <- all_paths %>% 
    filter(op == "~") %>% 
    mutate(from = lhs, to = rhs, style = "solid") %>% 
    select(from, to, style)
  edges_df <- combine_edges(loadings, regressions, predictions)
  
  graph <- create_graph(
    nodes_df = nodes_df,
    edges_df = edges_df,
    graph_attrs = c("layout = 'dot'")
  )
  return(graph)
}

先のlavaanで実行したfitfit2を実際に実行すると、以下のとおりです:

DgmR1 <- create_graph_lavaan(fit)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
## 
## Attaching package: 'Hmisc'
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      format.pval, round.POSIXt, trunc.POSIXt, units
## 
## Attaching package: 'dplyr'
##  以下のオブジェクトは 'package:Hmisc' からマスクされています: 
## 
##      combine, src, summarize
##  以下のオブジェクトは 'package:stats' からマスクされています: 
## 
##      filter, lag
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      intersect, setdiff, setequal, union
DgmR2 <- create_graph_lavaan(fit2)

render_graph(DgmR1)

render_graph(DgmR2)

課題

現在の状態では配置に問題があるため、改善が必要となります。またCFAと普通のSEMにしか対応していないので、他母集団SEMに対しても対応させたい。

以上です。